%Qt detection algotithm implemented from
%"New algorithm for QT interval analysis in 24-hour Holter ECG:
% perfrormance and applications"- P Laguna N.V. Thakor 
% @param x q dimension array with the ECG to analyze
% @param fs sampling frecuency of the ECG
% @return vector with the processed qt and qtp intervals
% @author Alex Mantaut
function [retQtIntervals, retQtpIntervals, qrs_on,t_off] = qt_detector(x, fs)

    close all;
    [differentiated, preprocessed] = qt_detector_filter(x, fs);
 
    pk1FindIntervalTime = 2;
    pk1FindIntervalTimeSamples = floor(pk1FindIntervalTime*fs);
    
    pknFindIntervalTime = 0.3;
    pknFindIntervalTimeSamples = floor(pknFindIntervalTime*fs);
    
    qpFindIntervalTime = 0.08;
    qpFindIntervalTimeSamples = floor(qpFindIntervalTime*fs);
    
    peakPrecedingFindIntervalTime = 0.1;
    peakPrecedingFindIntervalTimeSamples = floor(peakPrecedingFindIntervalTime*fs);
    
    % Peak before/after search times obtained from 
    % An algorithm for QRS onset and offset detection in single lead
    % electrocardiogram records
    pKaFindIntervalTime = 0.15;
    pKaFindIntervalTimeSamples = floor(pKaFindIntervalTime*fs);
    
    pKbFindIntervalTime = 0.3;
    pKbFindIntervalTimeSamples = floor(pKbFindIntervalTime*fs);
    
    rWaveFindIntervalTime = 0.3;
    rWaveFindIntervalTimeSamples = floor(rWaveFindIntervalTime*fs);
    
    maxQrDistance = 0.08;
    maxQrDistanceSamples = floor(maxQrDistance*fs);
    
    qrs1FindIntervalTime = 0.2;
    qrs1FindIntervalTimeSamples = floor(qrs1FindIntervalTime*fs);
    
    peakInhibitSearchTime = 0.15;
    peakInhibitSearchTimeSamples = floor(peakInhibitSearchTime*fs);
    
    %Time for learning of the qt detector (in qt intervals)
    qtLearningInterval = 6;
    qtpLearningInterval = 6;
    
    pKArray = zeros(1, length(preprocessed));
    RRArray = zeros(1, length(pKArray));
    qWaves = zeros(1, length(pKArray));
    rWaves = zeros(1, length(pKArray));
    qrs_on = zeros(1, length(pKArray));
    t_off = zeros(1, length(pKArray));
    qtIntervalsTime = zeros(1, length(pKArray));
    qtpIntervalsTime = zeros(1, length(pKArray));
    
    Kq = 2;
    Kr = 5;
    
    t = [0 :(1/fs) :(1/fs)*(length(rWaves)-1)];
   

    % Search for pk1
    [pK1AnalizedInterval,intervalStart,intervalEnd]  = subinterval(preprocessed,1,pk1FindIntervalTimeSamples);
    [pK1Value,pK1RelativePosition] =  max(abs(pK1AnalizedInterval));
    RRav = 0.8;
%     intervalStart =0;
    
    pKnPosition = 1;
    Hpeak = 0.8 * pK1Value;
    i=1;
    qtIntervalIterator = 1;
    qtpIntervalIterator = 1;
    
%     out = preprocessed;
    
    %---------------------detect R waves-----------------------------
    while (i <length(preprocessed))
        
        peakDetected = 0;
        RRn = (i - pKnPosition)/fs;
        
        if(abs(preprocessed(i)) >Hpeak)
            peakDetected = 1;
            pKnprevPosition = pKnPosition;
           
            [pKnAnalizedInterval,intervalStart,intervalEnd] = subinterval(abs(preprocessed),i,i+pknFindIntervalTimeSamples);
            
%             qWaves(intervalStart) = 0.3;
%             qWaves(intervalEnd) = 0.3;
            
            [peakValue,peakPosition] = max(pKnAnalizedInterval);
            
            pKnPosition = intervalStart +peakPosition -1;
        else
            %if no peak detected, go backwards
            if( RRn > 1.8*RRav)
                Hpeak = 0.8 * Hpeak;
                i = pKnPosition+ peakInhibitSearchTimeSamples;
                searchingBackwards = 1;
            end
        end
        
        if(peakDetected == 1)
            
            if(RRn < 0.6)
                peakDetected = 1;
            end
            
            searchingBackwards = 0;
            if( RRn < 1.5*RRav & RRn > 0.5*RRav)
               RRav = 0.8 * RRav +0.2 * RRn;
            end
            
            Hpeak = 0.8 *Hpeak + 0.2*abs(preprocessed(pKnPosition)); 
            pKArray(pKnPosition) = 1;
            RRArray(pKnPosition) = RRn;
            
            if( i +peakInhibitSearchTimeSamples < length(preprocessed))
                i = i + peakInhibitSearchTimeSamples;
            else
                i = length(preprocessed);
            end
        end
        
        i = i+1;
    end
    %----------------------------------------------------------------
    
    %todo- ver si puedo unificar los ciclos
    pKnPosition = 1;
%     RRav = 0.8;
    
    for i=1:length(pKArray)
         % if found a r wave
         if(pKArray(i) > 0)
             
             RRn = (i - pKnPosition)/fs;
             if( RRn < 1.5*RRav & RRn > 0.5*RRav)
               RRav = 0.8 * RRav +0.2 * RRn;
             end
             
             pKnPosition = i;
             
             %--------------------Rwave detection--------------------------
             [pKaAnalizedInterval,pKaIntervalStart, pKaIntervalEnd] = subinterval(preprocessed,pKnPosition,pKnPosition+pKaFindIntervalTimeSamples);
             pKaExtremes =  abspeakdet(pKaAnalizedInterval, 0.00005);
             
             peakDetectedSize = size(pKaExtremes);

             if peakDetectedSize(1) > 0
                 %In case there are no more peaks, It should take the pkn
                 if(peakDetectedSize(1) > 1)
                    pKaRelativePosition = pKaExtremes(2,1);
                 else
                    pKaRelativePosition = pKaExtremes(1,1);
                 end
             else
                 %HACK in case of arriving to the end of the array 
                 %pKaRelativePosition = 1;
                 break;
             end   
                 
             pKaValue = pKaAnalizedInterval(pKaRelativePosition);
              
             [pKbAnalizedInterval,pKbIntervalStart, pKbIntervalEnd] = subinterval(preprocessed,pKnPosition-pKbFindIntervalTimeSamples,pKnPosition);
             pKbExtremes =  abspeakdet(pKbAnalizedInterval, 0.00005);
             pKbRelativePosition = pKbExtremes(end,1);
             pKbValue = pKbAnalizedInterval(pKbRelativePosition);
             
             %Rp= zero-crossing preceding pKn
             if(abs(pKbValue) > abs(pKaValue))
                 [rpAnalizedInterval,intervalStart,intervalEnd] = subinterval(preprocessed,pKnPosition-rWaveFindIntervalTimeSamples,pKnPosition);
                 zeroPositionArray = crossing(rpAnalizedInterval,[],0,'none');
                 
                if length(zeroPositionArray) > 0
                    rpRelativePosition = zeroPositionArray(end);
                else
                    [minValue,minPosition] = min(abs(rpAnalizedInterval));
                    rpRelativePosition = minPosition;
                end
                
                rpPosition = intervalStart+rpRelativePosition-1;
                
             %Rp= zero-crossing following pKn
             else              
                [rpAnalizedInterval,intervalStart,intervalEnd] = subinterval(preprocessed,pKnPosition,pKnPosition+rWaveFindIntervalTimeSamples);
                zeroPositionArray = crossing(rpAnalizedInterval,[],0,'none');
                
                if length(zeroPositionArray) > 0
                    rpRelativePosition = zeroPositionArray(1);
                else
                    [minValue,minPosition] = min(abs(rpAnalizedInterval));
                    rpRelativePosition = minPosition;
                end
                
                rpPosition = intervalStart+rpRelativePosition-1;

             end
             
             rWaves(pKnPosition) = 0.3;
             rWaves(rpPosition) = 0.2;
             
             %--------------------Qp detection-------------------------
             [qpAnalizedInterval,intervalStart,intervalEnd] = subinterval(differentiated,rpPosition-qpFindIntervalTimeSamples,rpPosition-1);

             zeroPositionArray = crossing(qpAnalizedInterval,[],0,'none');
             
             extremePositions = peakdet(qpAnalizedInterval,0.005);

%              if(isempty(extremePositions) == 0)
%                  pknPosition = intervalStart -1 + extremePositions(end,1);
% %                  qWaves(pknPosition) = 0.3;
%              end
             
             %If no zeros found, the Q wave was detected as an R-wave 
             if(isempty(zeroPositionArray))
                 qpPosition = rpPosition;
                 qWavePresent = 0;
             else
                 firstQpZero = intervalStart +zeroPositionArray(end) -1;
                 
                 %If distance between qp and rp exceeds limits,
                 %then no q wave present
                 if(rpPosition-firstQpZero > maxQrDistanceSamples)
                     qpPosition = rpPosition;
                 else
                     qpPosition = firstQpZero;
                 end
                 
                 qWavePresent = 1;
             end
             
%              qWaves(qpPosition) = 0.2;

             %--------------------Qi detection----------
             [peakPrecedingAnalizedInterval,intervalStart,intervalEnd] = subinterval(differentiated, qpPosition- peakPrecedingFindIntervalTimeSamples,qpPosition);
             
             peakPrecedingRelativePosition =  abspeakdet(peakPrecedingAnalizedInterval,0.000005);
               
             if(isempty(peakPrecedingRelativePosition))
                [peakPrecedingValue,peakPrecedingRelativePosition] =  max(abs(peakPrecedingAnalizedInterval));
                peakPrecedingPosition = intervalStart -1 + peakPrecedingRelativePosition;
             else
                peakPrecedingPosition = intervalStart -1 + peakPrecedingRelativePosition(end,1);
             end
             
%              qWaves(peakPrecedingPosition) = 0.1;

             %Define QRS onset detection threshold
             if(qWavePresent == 1)
                 
                 qiPosition = peakPrecedingPosition;
                 
                 %differentiated(qiPosition)
                 %differentiated(i)/10
                 
                 if(abs(differentiated(qiPosition)) > abs((differentiated(i)/10)))
                     Hthres = differentiated(qiPosition) / Kq;
                 else
                     qWavePresent = 0;
                 end
             end
             
             if(qWavePresent == 0)
                 riPosition = peakPrecedingPosition;
                 Hthres = differentiated(riPosition) / Kr;
             end
             
            
             
             %Find QRS1
             [qrs1AnalizedInterval,intervalStart,intervalEnd] = subinterval(differentiated,peakPrecedingPosition- qrs1FindIntervalTimeSamples,peakPrecedingPosition);
             %TODO que hago si el cruce es negativo
             qrs1CrossingArray =  crossing((qrs1AnalizedInterval),[],Hthres,'none');
             
             if(~isempty(qrs1CrossingArray))
                qrs1Position = intervalStart + qrs1CrossingArray(end)-1;
                %qWaves(qrs1Position) = Hthres;
                qWaves(qrs1Position) = 0.2;
                qrs_on(qrs1Position) = 1;
             else
                 [minValue,minPosition] = min(abs(qrs1AnalizedInterval));
                 qrs1RelativePosition = minPosition;
                 qrs1Position = intervalStart + qrs1RelativePosition-1;
             end
             
             %----------------Bwind/Ewind--------------------
             
             if(RRav >0.7)
                 bwind = 0.14;
                 ewind = 0.5;
             else
                 bwind = 0.1;
                 ewind = 0.7*RRav;                 
             end
             
             bwindSamples = floor(bwind*fs);
             ewindSamples = floor(ewind*fs);
             
             [analizedWindow,analizedWindowStart,analizedWindowEnd] = subinterval(preprocessed,rpPosition+bwindSamples,rpPosition+ewindSamples);
             [windowMaxValue,windowMaxPosition] = max(analizedWindow);
             [windowMinValue,windowMinPosition] = min(analizedWindow);
             
%              qWaves(analizedWindowStart) = 0.5;
%              qWaves(analizedWindowEnd) = 0.5;
             rWaves(intervalStart+windowMaxPosition) = windowMaxValue;
             rWaves(intervalStart+windowMinPosition) = windowMinValue;
             
             
             %Create arrays with all the extremes
             windowExtremeValueArray = [windowMaxValue windowMinValue];
             windowExtremePositionArray = [windowMaxPosition windowMinPosition];
             windowExtremeIntervalStart = [analizedWindowStart analizedWindowStart];

             %-------------T wave Morphology detection----------
             % type 1: normal T wave (up-down)
             % type 2: inverted T wave (down-up)
             % type 3: only down T wave (down-down)
             % type 4: only up T wave (up-up)
             
             
             if(windowMaxPosition < windowMinPosition)
                 %type 1 or type 4
                 if(abs(windowMaxValue) < 4 *abs(windowMinValue))
                     tWaveType = 1;
                 else
                     tWaveType = 4;
                 end
             else
                 [minaAnalizedWindow,minaIntervalStart,minaIntervalEnd] = subinterval(analizedWindow,windowMaxPosition,windowMaxPosition+length(analizedWindow));
                 
                 [windowMinaValue,windowMinaPosition] = min(minaAnalizedWindow);
                 
                 %Add mina to the extreme arrays
                 windowExtremeValueArray(3) = windowMinaValue;
                 windowExtremePositionArray(3) = windowMinaPosition;
                 windowExtremeIntervalStart(3) = analizedWindowStart;
                 
                 if(abs(windowMaxValue) < 4 *abs(windowMinaValue))
                     tWaveType = 1;
                 else
                     if(abs(windowMinValue) > 4 *abs(windowMaxValue))
                         tWaveType = 3;
                     else
                         tWaveType = 2;
                     end
                 end
             end
             
             %search for Ti
             %look for the last extreme
             [tiLastExtremePosition,tiLastExtremeType] = max(windowExtremePositionArray);
             tiPosition = windowExtremeIntervalStart(tiLastExtremeType) -1 + tiLastExtremePosition;
             tiValue = abs(preprocessed(tiPosition));
             qWaves(tiPosition) = 0.2;
             
             %Define theshold for T2
             Ht = tiValue/2;
             
             [t2AnalizedInterval,t2IntervalStart,intervalEnd] = subinterval(preprocessed,tiPosition,analizedWindowEnd);
             t2CrossingArray =  crossing(abs(t2AnalizedInterval),[],Ht,'none');
             
             
             [t1AnalizedInterval,t1IntervalStart,intervalEnd] = subinterval(preprocessed,analizedWindowStart,tiPosition);
             t1CrossingArray =  crossing(t1AnalizedInterval,[],0,'none');
             
             t1Position = 0;
             t2Position = 0;
             
             
             %QT interval calculation
             if(~isempty(t2CrossingArray))
                t2Position = t2IntervalStart+t2CrossingArray(1)-1;
                qWaves(t2Position) = 0.2;
                t_off(t2Position) = 1;
                
                currentQtInterval = (t2Position - qrs1Position)/(fs*sqrt(RRn));
               
                qtIntervals(qtpIntervalIterator) = currentQtInterval;  
                qtIntervalIterator = qtIntervalIterator +1;
                qtIntervalsTime(i) = currentQtInterval;
             end
             

             %calculate qtp
            if(~isempty(t1CrossingArray))
                t1Position = t1IntervalStart+t1CrossingArray(end)-1;
%                 qWaves(t1Position) = 0.15;

                currentQtpInterval = (t1Position - qrs1Position)/(fs*sqrt(RRn));
                
                qtpIntervals(qtpIntervalIterator) = currentQtpInterval;  
                qtpIntervalIterator = qtpIntervalIterator +1;
                qtpIntervalsTime(i) = currentQtpInterval;
            end
         end     
    end
    
    retQtpIntervals = qt_detector_post_processing(qtpIntervals);
    retQtIntervals = qt_detector_post_processing(qtIntervals);
    
    senAux = (x - mean(x))/max(x);
    senAux2 = differentiated; % /max(differentiated);
    senAux3 = preprocessed; %/max(preprocessed);
    
    rWaves= subinterval(rWaves, 1, length(pKArray));
    
%     figure(1);
%     hold;
%     signalPlot = plot(t,senAux,'r');
%     signalPlot = plot(t,senAux3,'c');
%     signalPlot = plot(t,senAux2,'g');
%     rPlot = stem(t,rWaves,'.r');
    %legend(rPlot,'R waves');
%     qPlot = stem(t,qWaves,'.b');
%     legend(qPlot,'Q waves');

%     figure(2);
%     hold;
%     plot(t,senAux,'r');
%     plot(t,senAux2,'g');
% %     plot(t,senAux3,'c');
%     title('Detected points')
%     qPlot = stem(t,qWaves,'.b');
% %     qPlot = stem(t,rWaves,'.r');
% %     legend(qPlot,'Q waves');
    
%         figure(2);
%     hold;
%     plot(senAux,'r');
%     plot(senAux2,'g');
% %     plot(t,senAux3,'c');
%     title('Detected points')
%     stem(qWaves,'.b');
% %     stem(rWaves,'.r');
%     stem(pKArray/5,'.y');
%     
%    
%     figure(3)
%     subplot(4,1,1);
%     hist(qtIntervals,40);   
%     title('QT Intervals')
%     subplot(4,1,2);
%     hist(retQtIntervals,40);
%     subplot(4,1,3);
%     stem(qtIntervalsTime);
%     subplot(4,1,4);
%     stem(RRArray);
%     
%     figure(4)
%     subplot(3,1,1);
%     hist(qtpIntervals,40);
%     title('QTp Intervals')
%     subplot(3,1,2);
%     hist(retQtpIntervals,40);
%     subplot(3,1,3);
%     stem(qtpIntervalsTime);

end

function retIntervals = qt_detector_post_processing(inputIntervals) 

%eliminate the first and last meassurements (possible incomplete beats) 
analizedInterval= subinterval(inputIntervals,2,length(inputIntervals)-1);
retIntervals =  analizedInterval;

learningInterval = 5;
toBeRemoved = [];

for i=1:length(analizedInterval)
    %if the learning interval has passed
    if i > learningInterval
        [analizedSubInterval,intervalStart,intervalEnd] = subinterval(analizedInterval,1,i);
        intervalMean =  mean(analizedSubInterval);     
 
        %Check if the new beat is less than 15% than the
        %average qt
        if abs(analizedInterval(i) -intervalMean) > 0.15 *  intervalMean
            toBeRemoved(end+1) = i;
        end
    end
end

retIntervals(toBeRemoved) = [];

analizedInterval = retIntervals;

toBeRemoved = [];

for i=1:5:length(analizedInterval)
    [analizedSubInterval,intervalStart,intervalEnd] = subinterval(analizedInterval,i,i+5);
    
    [maxValue,maxPosition] = max(analizedSubInterval);
    [minValue,minPosition] = min(analizedSubInterval);
    toBeRemoved(end+1) = intervalStart + minPosition -1;
    toBeRemoved(end+1) = intervalStart + maxPosition -1;
end

retIntervals(toBeRemoved) = [];

end
